// compiles as:
//
// gcc fftw_guru_test.c -L/usr/local/fftw/fftw-3.3.3/lib  -lfftw3 -lm 
// -I/usr/local/fftw/fftw-3.3.3/include -o c_test.x

#include <fftw3.h>   

int main()
{
  const int nx = 8;               // 3D array dimenrions    
  const int ny = 6;
  const int nz = 4;

  const int jx = 3;               // point at which 1D and 3D 
  const int jy = 3;               // transforms to be compared
  const int jz = 3;


  fftw_complex   u[nz][ny][nx];   // input 3D array
  fftw_complex   v[nz][ny][nx];   // output 3D array

  fftw_complex   ux[nx];          // in-place 1D arrays 
  fftw_complex   uy[ny];          // for comparison
  fftw_complex   uz[nz];


  int ix, iy, iz;                 // running indices

  fftw_plan p1x, p1y, p1z, p3x, p3y, p3z;

  int                rank = 1;
  fftw_iodim         dims[1];
  int                hm_rank = 2; 
  fftw_iodim         hm_dims[2];
 

  /*-- initiaization of data --*/

  for (iz = 0; iz<nz; iz++) 
    for (iy = 0; iy<ny; iy++) 
      for (ix = 0; ix<nx; ix++) {
	u[iz][iy][ix][0] = drand48();
        u[iz][iy][ix][1] = 0;
      }


  /*--  1D transforms of pensils passing through jx,jy,iz point --*/

  for (ix = 0; ix<nx; ix++){
    ux[ix][0] = u[jz][jy][ix][0];
    ux[ix][1] = u[jz][jy][ix][1];
  }

  for (iy = 0; iy<ny; iy++){
    uy[iy][0] = u[jz][iy][jx][0];
    uy[iy][1] = u[jz][iy][jx][1];
  }

  for (iz = 0; iz<nz; iz++){
    uz[iz][0] = u[iz][jy][jx][0];
    uz[iz][1] = u[iz][jy][jx][1];
  }

  p1x = fftw_plan_dft_1d(nx, ux, ux, FFTW_FORWARD, FFTW_ESTIMATE);
  p1y = fftw_plan_dft_1d(ny, uy, uy, FFTW_FORWARD, FFTW_ESTIMATE);
  p1z = fftw_plan_dft_1d(nz, uz, uz, FFTW_FORWARD, FFTW_ESTIMATE);
 
  fftw_execute(p1x);
  fftw_execute(p1y);
  fftw_execute(p1z);


 /*--  1D transforms of 3D array in x-direction --*/

  dims[0].n     = nx;
  dims[0].is    = dims[0].os = 1;

  hm_dims[1].n  = ny;
  hm_dims[1].is = hm_dims[1].os = nx;

  hm_dims[0].n  = nz;
  hm_dims[0].is = hm_dims[0].os = nx*ny;

 
  p3x = fftw_plan_guru_dft(rank, dims, hm_rank, hm_dims,
			   &u[0][0][0], &v[0][0][0],  
			   FFTW_FORWARD, FFTW_ESTIMATE);

  fftw_execute(p3x);

  printf("\n Testing guru interface in x-direction:\n\n");

  for (ix = 0; ix < nx; ix++) 
    printf("%10.3e + i ( %10.3e )       %10.3e + i ( %10.3e )\n",
	   ux[ix][0],  ux[ix][1],   v[jz][jy][ix][0],  v[jz][jy][ix][1] );  



 /*--  1D transforms of 3D array in y-direction --*/

  dims[0].n     = ny;
  dims[0].is    = dims[0].os = nx;

  hm_dims[1].n  = nx;
  hm_dims[1].is = hm_dims[1].os = 1;

  hm_dims[0].n  = nz;
  hm_dims[0].is = hm_dims[0].os = nx*ny;
 

  p3y = fftw_plan_guru_dft(rank, dims, hm_rank, hm_dims,
			   &u[0][0][0], &v[0][0][0],  
			   FFTW_FORWARD, FFTW_ESTIMATE);

  fftw_execute(p3y);

  printf("\n\n Testing guru interface in y-direction:\n\n");

  for (iy = 0; iy < ny; iy++) 
    printf("%10.3e + i ( %10.3e )       %10.3e + i ( %10.3e )\n",
	   uy[iy][0],  uy[iy][1],   v[jz][iy][jx][0],  v[jz][iy][jx][1] );  



 /*--  1D transforms of 3D array in z-direction --*/

  dims[0].n     = nz;
  dims[0].is    = dims[0].os = nx*ny;

  hm_dims[1].n  = nx;
  hm_dims[1].is = hm_dims[1].os = 1;

  hm_dims[0].n  = ny;
  hm_dims[0].is = hm_dims[0].os = nx;
 

  p3z = fftw_plan_guru_dft(rank, dims, hm_rank, hm_dims,
			   &u[0][0][0], &v[0][0][0],  
			   FFTW_FORWARD, FFTW_ESTIMATE);

  fftw_execute(p3z);

  printf("\n\n Testing guru interface in z-direction:\n\n");

  for (iz = 0; iz < nz; iz++) 
    printf("%10.3e + i ( %10.3e )       %10.3e + i ( %10.3e )\n",
	   uz[iz][0],  uz[iz][1],   v[iz][jy][jx][0],  v[iz][jy][jx][1] );  


  /*-- done --*/

  printf("\n\n Two columns of complex numbers should match.\n\n");

  return(0);

}

